Starting to look at linear regression for predicting pitch pine presence/abundance. I think markdown scripts are easier to keep neat and organize.
All counts are included in this script (PIRI, Shrub Oak, and Other categories)
Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Functions:
# functions -------------------
Import small seedling data:
# global variables -------------------
source("Scripts/SS_Import.R")
Pivot wider to create dataframe where each row is for one plot and has total details for each species group
ss_merge2 <- ss_merge %>%
select(Plot_No, Region, Treat_Type, Site, Species_Groups, Total) #this is dropping stump sprout, browse, and germinate data
ss_merge2 <- ss_merge2 %>%
pivot_wider(names_from = Species_Groups, values_from = Total)
Import time since data and add it to the small seedling dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
ss_merge3 <- merge(ss_merge2, time_since, by = "Site")
#log transform time from treatment data
ss_merge3$l.TFT <- log(ss_merge3$Time_from_Treat)
Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with ss dataset -------------------
ss_merge4 <- merge(ss_merge3, prism_BA, by = "Plot_No")
Run ‘Ground_Data.R’ script and add it to small seedling dataset:
source("Scripts/Ground_Data.R")
# merge with ss dataset -------------------
ss_merge5 <- merge(ss_merge4, ground3, by = "Plot_No")
rm(ss_merge2,
ss_merge3,
ss_merge4,
time_since)
Source and add veg data
source("Scripts/Veg_Data.R")
# merge with ss dataset
ss.all <- merge(ss_merge5, veg3, by = "Plot_No")
Create log transformed categories for newly added variables, then select for just the desired variables:
ss.all$l.PIRI <- log(ss.all$PIRI + 1)
ss.all$l.SO <- log(ss.all$Shrub_Oak + 1)
ss.all$l.other <- log(ss.all$Other + 1)
ss.all2 <- ss.all %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)
Select just for numerical vs log and then look at paired plots:
#not transformed
ss.num <- ss.all2 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(ss.num, aes(color = Treat_Type))
ggpairs(ss.num)
#log transformed
ss.numl <- ss.all2 %>%
select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(ss.numl, aes(color = Treat_Type))
ggpairs(ss.numl)
rm(ss.num,
ss.numl)
Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)
Full dataset is called ss.all2.
Hidden below is the original analysis and step by step model elimination conducted without treatment type included. Best model without treatment type had a higher AIC and is underdispersed according to DHARMa tests
Model step wise elimination with treatment type included
ss.m1 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m1) #724.86
## [1] 724.8561
#would like to test ba vs ba piri
ss.m2a <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m2a) #723.6
## [1] 723.6111
lrtest(ss.m1, ss.m2a) # p = .4
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -348.43
## 2 13 -348.81 -1 0.7551 0.3849
ss.m2b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m2b) #725.9
## [1] 725.919
lrtest(ss.m1, ss.m2b) # p =.08, so not important
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -348.43
## 2 13 -349.96 -1 3.0629 0.0801 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m2a, ss.m2b)
# test ld vs mineral exposure
ss.m3a <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m3a) #723.0
## [1] 722.9976
lrtest(ss.m1, ss.m3a) #p=.25
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -348.43
## 2 11 -350.50 -3 4.1415 0.2466
ss.m3b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m3b) #730.1
## [1] 730.1247
lrtest(ss.m1, ss.m3b) # p = 0.01, loss of avg ld significant
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -348.43
## 2 11 -354.06 -3 11.269 0.01036 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m1, ss.m3b)
# test l.other
ss.m4 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m4) #722.6
## [1] 722.5895
lrtest(ss.m3a, ss.m4) #p=.2, can drop other
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -350.50
## 2 10 -351.29 -1 1.592 0.207
rm(ss.m3a)
# test l.SO
ss.m5 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m5) #725.5
## [1] 725.4735
lrtest(ss.m4, ss.m5) #p = 0.03, l.so important
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -351.29
## 2 9 -353.74 -1 4.884 0.02711 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m5)
#test veg
ss.m6 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(ss.m6) #727
## [1] 726.9952
lrtest(ss.m4, ss.m6) # p = 0.011, veg important
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -351.29
## 2 9 -354.50 -1 6.4056 0.01138 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(ss.m6)
do pairwise on variables of interest from reduced model
ss.pair <- ss.all2 %>%
select(Region, Treat_Type, l.PIRI, l.SO, l.Veg_Total, avgLD_l, l.TFT)
ggpairs(ss.pair, aes(color = Treat_Type))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rm(ss.pair)
Not seeing any clear, strong correlations
Best model:
ss.m4 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
summary(ss.m4) #722.6
## Family: poisson ( log )
## Formula:
## PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Data: ss.all2
##
## AIC BIC logLik deviance df.resid
## 722.6 764.7 -351.3 702.6 488
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 4.657 2.158
## Site (Intercept) 2.247 1.499
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.2307 1.7872 -2.927 0.00343 **
## Treat_TypeFallRx 6.8811 1.4362 4.791 0.00000166 ***
## Treat_TypeHarvest 6.4777 1.5012 4.315 0.00001595 ***
## Treat_TypeMowRx 5.2753 1.4183 3.719 0.00020 ***
## Treat_TypeSpringRx 6.7837 1.4407 4.709 0.00000249 ***
## l.SO -0.5172 0.2390 -2.164 0.03049 *
## avgLD_l -1.1967 0.4340 -2.757 0.00583 **
## l.Veg_Total -0.8041 0.3280 -2.451 0.01423 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#This is exactly the same, just checking Maria's comment
ss.m4a <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site) + (1|Plot_No),
data = ss.all2,
family = poisson)
summary(ss.m4a)
## Family: poisson ( log )
## Formula:
## PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site) + (1 | Plot_No)
## Data: ss.all2
##
## AIC BIC logLik deviance df.resid
## 722.6 764.7 -351.3 702.6 488
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 2.247 1.499
## Plot_No (Intercept) 4.657 2.158
## Number of obs: 498, groups: Site, 47; Plot_No, 498
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.2307 1.7872 -2.927 0.00343 **
## Treat_TypeFallRx 6.8811 1.4362 4.791 0.00000166 ***
## Treat_TypeHarvest 6.4777 1.5012 4.315 0.00001595 ***
## Treat_TypeMowRx 5.2753 1.4183 3.719 0.00020 ***
## Treat_TypeSpringRx 6.7837 1.4407 4.709 0.00000249 ***
## l.SO -0.5172 0.2390 -2.164 0.03049 *
## avgLD_l -1.1967 0.4340 -2.757 0.00583 **
## l.Veg_Total -0.8041 0.3280 -2.451 0.01423 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Checking model fit:
ss.m4_sr <- simulateResiduals(ss.m4, n = 1000, plot = TRUE) #looks pretty good
testResiduals(ss.m4_sr) #passes uniformity, dispersion, and outliers
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055502, p-value = 0.09301
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0070259, p-value = 0.212
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00118473895582329 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.055502, p-value = 0.09301
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0070259, p-value = 0.212
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00118473895582329 )
## 0
testQuantiles(ss.m4_sr) #passes
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.2776
## alternative hypothesis: both
testZeroInflation(ss.m4_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99615, p-value = 0.946
## alternative hypothesis: two.sided
Time to pairwise test
emmeans(ss.m4, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type rate SE df asymp.LCL asymp.UCL
## Control 0.000147 0.000196 Inf 0.0000108 0.002
## FallRx 0.143092 0.098074 Inf 0.0373430 0.548
## Harvest 0.095594 0.081151 Inf 0.0181066 0.505
## MowRx 0.028723 0.019958 Inf 0.0073586 0.112
## SpringRx 0.129822 0.093714 Inf 0.0315422 0.534
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0.00103 0.00147 Inf 1 -4.791 <.0001
## Control / Harvest 0.00154 0.00231 Inf 1 -4.315 0.0002
## Control / MowRx 0.00512 0.00726 Inf 1 -3.719 0.0019
## Control / SpringRx 0.00113 0.00163 Inf 1 -4.709 <.0001
## FallRx / Harvest 1.49687 1.48518 Inf 1 0.407 0.9943
## FallRx / MowRx 4.98172 4.23482 Inf 1 1.889 0.3232
## FallRx / SpringRx 1.10222 0.97175 Inf 1 0.110 1.0000
## Harvest / MowRx 3.32809 3.28617 Inf 1 1.218 0.7411
## Harvest / SpringRx 0.73635 0.73882 Inf 1 -0.305 0.9981
## MowRx / SpringRx 0.22125 0.18896 Inf 1 -1.766 0.3935
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
FallRx, Harvest, MowRx, SpringRx all significantly different from Control - no treatment significantly different from another
Model ss.m4 uses Treatment Type, Scrub Oak presence, average leaf litter depth, and total veg cover as predictors. Make a graph for each
I’d like to see some graphs before I close out
library(sjPlot)
library(sjlabelled)
##
## Attaching package: 'sjlabelled'
## The following object is masked from 'package:forcats':
##
## as_factor
## The following object is masked from 'package:ggplot2':
##
## as_label
## The following object is masked from 'package:dplyr':
##
## as_label
library(sjmisc)
##
## Attaching package: 'sjmisc'
## The following object is masked from 'package:purrr':
##
## is_empty
## The following object is masked from 'package:tidyr':
##
## replace_na
## The following object is masked from 'package:tibble':
##
## add_case
set_theme(base = theme_classic(),
theme.font = 'serif',
axis.title.size = 1.5,
axis.textsize.x = 1.5,
axis.textsize.y = 1.5,
title.size = 2.5,
title.align = "center",
legend.pos = "right",
legend.size = 1.5,
legend.title.size = 1.5,
#legend.bordercol = "black",
legend.item.size = .75)
plot_model(ss.m4) #this is incidence rate ratios
plot_model(ss.m4, type = "diag") #random vs normal quantiles
## $`Plot_No:Site`
## `geom_smooth()` using formula = 'y ~ x'
##
## $Site
## `geom_smooth()` using formula = 'y ~ x'
plot_model(ss.m4, type = "re") #this plots random effects
## [[1]]
##
## [[2]]
plot_model(ss.m4,
type = 'pred',
terms = 'Treat_Type') #plot marginal effects (i might have done this wrong)
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
# SS PIRI plot 1 ------------------- LD vs. TT
plot_model(ss.m4,
type = 'pred',
terms = c('avgLD_l', 'Treat_Type'),
axis.title = c("Average Leaf Litter Depth (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Seedlings <50cm",
legend.title = "Treatment Type",
line.size = 1,
value.offset = 'Treat_Type',
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
# SS PIRI plot 2 ------------------- Veg vs TT
plot_model(ss.m4,
type = 'pred',
terms = c('l.Veg_Total', 'Treat_Type'),
axis.title = c("Average Understory Vegetation Cover (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Seedlings <50cm",
legend.title = "Treatment Type",
line.size = 1,
value.offset = 'Treat_Type',
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
# SS PIRI pot 3 ------------------- Shrub Oak seedling vs TT
plot_model(ss.m4,
type = 'pred',
terms = c('l.SO', 'Treat_Type'),
axis.title = c("Shrub oak small seedling counts (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Seedlings <50cm",
legend.title = "Treatment Type",
line.size = 1,
value.offset = 'Treat_Type',
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
# SS PIRI plot 4 ------------------- Time from Treatment vs TT
plot_model(ss.m4,
type = 'pred',
terms = c('l.TFT', 'Treat_Type'),
axis.title = c("Time from Treatment (log transformed)", "Total Count of Pitch Pine"),
title = "Predicted Counts of Pitch Pine Seedlings <50cm",
legend.title = "Treatment Type",
line.size = 1,
value.offset = 'Treat_Type',
ci.lvl = NA,
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'
Next steps:
2e. other interaction terms to think about? 3. take another look at treat type glmms 4. looks at lr sans treat type. Better AIC or fit? stronger predictors? maybe in LS or SA models? 6. start with full models post scatter plot analysis 7. are there other ways that I want to think about veg cover? e.g., just pulling out SO cover? VA/GA cover?
Begin model elimination
so.ss1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss1) #1897.4
## [1] 1897.374
#test piri BA
so.ss2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss2) #1896
## [1] 1895.956
# test total ba
so.ss3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss3) #1894.1
## [1] 1894.07
# test mineral
so.ss4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss4) #1892.3
## [1] 1892.281
# test litter depth
so.ss5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss5) #1890.1
## [1] 1890.831
# test other seedling
so.ss6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss6) #1889.5
## [1] 1889.516
lrtest(so.ss5, so.ss6) # p = .4
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -935.42
## 2 9 -935.76 -1 0.685 0.4079
# test PIRI seedling
so.ss7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss7) #1892.6
## [1] 1892.648
lrtest(so.ss6, so.ss7) # p = 0.02, so l.PIRI is important
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -935.76
## 2 8 -938.32 -1 5.1315 0.0235 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# add PIRI back in, test veg cover
so.ss8 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss8) #1888.9
## [1] 1888.943
lrtest(so.ss6, so.ss8) # p = 0.2
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -935.76
## 2 8 -936.47 -1 1.4263 0.2324
# test treat type
so.ss9 <- glmmTMB(Shrub_Oak ~ l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
AIC(so.ss9) #1911.5
## [1] 1911.528
lrtest(so.ss9, so.ss8) # p > 0.001, so TT is important
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -951.76
## 2 8 -936.47 4 30.585 0.000003719 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Best model, then test
so.ss8 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = poisson)
summary(so.ss8)
## Family: poisson ( log )
## Formula:
## Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ss.all2
##
## AIC BIC logLik deviance df.resid
## 1888.9 1922.6 -936.5 1872.9 490
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 1.923 1.387
## Site (Intercept) 1.485 1.219
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3474 0.4346 -10.003 < 0.0000000000000002 ***
## Treat_TypeFallRx 3.1111 0.6183 5.032 0.0000004851 ***
## Treat_TypeHarvest 1.3851 0.7377 1.877 0.0605 .
## Treat_TypeMowRx 3.1769 0.5759 5.516 0.0000000346 ***
## Treat_TypeSpringRx 2.9035 0.6316 4.597 0.0000042769 ***
## l.PIRI -0.4518 0.1890 -2.390 0.0168 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ss8_sr <- simulateResiduals(so.ss8, n = 1000, plot = TRUE) #passes KS, quantile deviations signficant
testResiduals(so.ss8_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033251, p-value = 0.6406
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0083534, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.96
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00184738955823293 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033251, p-value = 0.6406
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0083534, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.96
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00184738955823293 )
## 0
testDispersion(so.ss8_sr, alternative = "less") #fails, so it is underdispersed
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0083534, p-value = 0.001
## alternative hypothesis: less
#try again with gen pois / com pois ditros
so.ss10 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = compois)
summary(so.ss10)
## Family: compois ( log )
## Formula:
## Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ss.all2
##
## AIC BIC logLik deviance df.resid
## 1869.6 1907.5 -925.8 1851.6 489
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 1.191 1.091
## Site (Intercept) 1.504 1.226
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for compois family (): 4.39
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0292 0.4373 -9.213 < 0.0000000000000002 ***
## Treat_TypeFallRx 3.0398 0.6152 4.941 0.0000007773 ***
## Treat_TypeHarvest 1.4196 0.7350 1.931 0.0534 .
## Treat_TypeMowRx 3.1285 0.5737 5.453 0.0000000494 ***
## Treat_TypeSpringRx 2.9044 0.6288 4.619 0.0000038500 ***
## l.PIRI -0.4238 0.1779 -2.382 0.0172 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ss10_sr <- simulateResiduals(so.ss10, n = 1000, plot = TRUE)
testResiduals(so.ss10_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.040131, p-value = 0.3989
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.050676, p-value = 0.01
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00146586345381526 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.040131, p-value = 0.3989
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.050676, p-value = 0.01
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.00146586345381526 )
## 0
testDispersion(so.ss10_sr, alternative = "less") # still under dispersed
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.050676, p-value = 0.005
## alternative hypothesis: less
so.ss11 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ss.all2,
family = genpois)
summary(so.ss11)
## Family: genpois ( log )
## Formula:
## Shrub_Oak ~ Treat_Type + l.PIRI + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ss.all2
##
## AIC BIC logLik deviance df.resid
## 1854.2 1892.1 -918.1 1836.2 489
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.2882 0.5369
## Site (Intercept) 0.9904 0.9952
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for genpois family (): 5.87
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.4252 0.3608 -9.493 < 0.0000000000000002 ***
## Treat_TypeFallRx 2.8087 0.5032 5.581 0.000000023864 ***
## Treat_TypeHarvest 1.4256 0.6233 2.287 0.0222 *
## Treat_TypeMowRx 2.9454 0.4699 6.268 0.000000000366 ***
## Treat_TypeSpringRx 2.8173 0.5213 5.404 0.000000065125 ***
## l.PIRI -0.3657 0.1623 -2.253 0.0242 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.ss11_sr <- simulateResiduals(so.ss11, n = 1000, plot = TRUE)
testResiduals(so.ss11_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.025101, p-value = 0.9123
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.28065, p-value = 0.092
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00166666666666667 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.025101, p-value = 0.9123
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.28065, p-value = 0.092
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.008032129
## sample estimates:
## outlier frequency (expected: 0.00166666666666667 )
## 0
testDispersion(so.ss11_sr, alternative = "less") # at 0.05, so could possibly work
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.28065, p-value = 0.046
## alternative hypothesis: less
testQuantiles(so.ss11_sr) #quantiles are not good
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.00001043
## alternative hypothesis: both
emmeans(so.ss11, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.183 0.0666 Inf 0.0898 0.373
## FallRx 3.038 1.1048 Inf 1.4896 6.196
## Harvest 0.762 0.4032 Inf 0.2701 2.149
## MowRx 3.483 1.1174 Inf 1.8574 6.532
## SpringRx 3.064 1.2098 Inf 1.4134 6.643
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0.0603 0.0303 Inf 1 -5.581 <.0001
## Control / Harvest 0.2404 0.1498 Inf 1 -2.287 0.1489
## Control / MowRx 0.0526 0.0247 Inf 1 -6.268 <.0001
## Control / SpringRx 0.0598 0.0312 Inf 1 -5.404 <.0001
## FallRx / Harvest 3.9873 2.5162 Inf 1 2.192 0.1827
## FallRx / MowRx 0.8722 0.4141 Inf 1 -0.288 0.9985
## FallRx / SpringRx 0.9914 0.5212 Inf 1 -0.016 1.0000
## Harvest / MowRx 0.2188 0.1325 Inf 1 -2.509 0.0885
## Harvest / SpringRx 0.2486 0.1605 Inf 1 -2.156 0.1966
## MowRx / SpringRx 1.1367 0.5632 Inf 1 0.259 0.9990
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
FallRx, MowRx, and SpringRx significantly different than control; none of the treatments significantly different from each other
This analysis can be seen in the pair plots at the top of script
Veg cover, average litter depth, and time from treatment all have a clear negative relationship (as they decrease, total increases), some more significant than others
Basal area per hectare and mineral soil exposure have positive relationship (as they increase, total increase), these relationships appear weaker than the above negative relationships
More scatterplots hidden below; see pairwise at beginning of script